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Abstract: In this work, the equations of motion of the restricted three-body problem under the effects of the 
oblateness of less massive primary and the radiation pressure of the bigger massive primary are expressed. 
The analytical solution is obtained by the variation of parameters method. The locations of the libration 
points are obtained. The periodic orbits around each of these points are investigated for the Sun-Earth 
system, the zero-velocity curves, the phase spaces, and the Poincare surface sections are presented for one of 
the collinear libration points and one of the non-linear libration points. The obtained results are compared 
with previous works such as Kunitsyn!3] and Schuerman!!4l, and good agreements with these results are 


found. 


Introduction 

The variation of parameters method is one of the important methods of theoretical and 
applied mathematics. The applications of this method are concerned mainly with the 
determination of maximum and minimum of certain expressions involving unknown 
functions. For orbital dynamics this method was introduced by the Swiss-born 
mathematician Leonhard Euler (1707 - 1783) and completed by the Italian-French 
mathematician Joseph-Louis Lagrange (1736 - 1813) Abell [1]and Efroimsky[2] Lagrange 
had put a reduced two-body problem as an unperturbed solution and had assumed that all 
perturbations come from the gravitational force. In the 20th century Vallado[3], celestial 
mechanics began to consider interactions that depend on both positions and velocities 


(relativistic corrections, atmospheric drag, inertial forces). Therefore, the method of 


The Restricted Three-Body Problem. 
©2020 Science Asia 
1/14 


2 M.N. ISMAIL, A.H. IBRAHIM, SAHAR H. YOUNIS, GHADA F. MOHAMDIEN 


variation of parameters used by Lagrange was extended to the situation with velocity- 
dependent forces Moulton [4]. Geometrically, this method is a representation of an orbit as 
a set of points, each of which is contributed by a member of some chosen family of curves 


C(K), where K stands for a set of constants that number a particular -curve within the 
family curves Lovettl5! 
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Fig.l. Each point of the orbit is contributed by a member of some family of curves C(K) of a 
certain type, K standing for a set of constants that member a particular curve within the family. 
Motion from A to B is, first, due to the motion along the curve C(K) from 4, to A, and, second, 
due to the fact that during this motion the curve itself was evolving from C (K1) to C (K3). 


This situation is depicted within the family of curves Fig.1.Point A of the orbit coincides 
with some pointd,on a curve C (Ai). Point B of the orbit coincides with pointA,on some 
other curve C (K2) of the same family, etc. This way, orbital motion from A to B becomes a 
superposition of motion along Cxfrom/,toA,and a gradual distortion of the curve Cxfrom 
the shape C (1) to the shape C (K2). Normally the curves Cxare chosen to be ellipses or 
hyperbolae to be managed analyzed Oberti[6]. Many analytical theories depend on the 
central body and the perturbs force, therefore it is needed to study the problem with many 
perturbing forces Vallado[7]. 

Therefore, the method of variation of parameters is useful to study the motion of a 
massless particle in restricted three-body problem which considered as one of the most 
important objects in astro-dynamics. There are many studies which treated the restricted 


three-body problem under different perturbing forces Srivastava[8]. But most important 
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perturbing forces are the radiation pressure and the oblateness of the primaries 
Efroimsky [2], George [9], Douskos[10], Sharma[11], Simmons[12], Kunitsyn[13] and 
Schuerman[14]. 


In this work, the equations of motion of the restricted three-body problem in the 
classical form are presented under the effects of the radiation pressure of the more massive 
body and the oblateness of less massive body. The variation of parameters method is used 
to obtain the analytical solution. From this solution, the libration points are obtained and 


the stability around each point is studied. 


Equations of Motion 

Using a barycentric-synodic coordinate system (X, Y, Z) and dimensionless variables, the 
equations of motion of a test particle in the circular restricted three-body problem under 
the effects of the oblateness of the small primary and the radiation pressure of the bigger 


primary can be expressed as 


X — 2nY = Uy (1.1) 
Y + 2nX = Uy (1.2) 
Ž = U; (1.3) 


Since the above system is rotating around the Z axis by a constant angular velocity, then Z 


= constant, so the system of Equations (1) is reduced to the system 


Ž — 2nř = Uy (2.1) 
Y + 2nX = Uy (2.2) 
Where 

_ ny. yay, G-A- u uA 
U => (X quse n ad (3) 
Where 


r =y- +y’, a0 1-414 y* 


U, = n?x — p (1-B) _ nr _ 3A wow) (4) 
1 2 2 
a-.))0-8)y Huy 3Any 
Uy = my R$ i8 (5) 


Where n is the mean motion of the smaller primary and given by 
rè- 


2 
E EL Ibrahim[15], reand rpare the equatorial and polar radii of 





n?^-1- 2A, where A = 


the oblate body, and R is the distance between the two primaries, is considered as unity 
Ismail[16]. 
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Since the radiation, pressure force Fr and the gravitational force Fgare acting on the body in 


opposite directions, so that 





Fg -Fr = Fg (1-B) (6) 

Where 

Fg and Fr are the gravitational force and radiation force respectively, 

B =F In dimensionless coordinates zm + m2= 1, 4 = Tm is the mass ratio of the 
system. 


Method of Variation of Parameters 

To solve the system of Equations (2) we use the method of variation of parameters 
Chen[17] and Palais[18], which needed to reduce the above system to first order, so that let 
u = DX,Du = D?X, v = DY,Dv = D?Y. Then the system of Equations (2) becomes 


Du = 2nv + Uy (7.1) 
Dv = —2nu + Uy (7.2) 
DX =u (7.3) 
DY =v (7.4) 
Equations (7) are system of ODEs, which expressed in the matrix form as 
ù 0 2n 0 07 pu Uy 
zoslal- Fe $8 8 jeu o 
Y 0 1 oo lY 0 
Now, the homogenous and the particular- solutions for the system of Equations (8) well be 
obtained. At first the homogenous solution yy is obtained when Ux= 0 and Uy- 
0;then 
ü 0 2n 0 Q fu 0 
E MEME: o 
Y 0 1 0041Y 0 
The auxiliary equation for the homogenous system of Equations (9) is 
0—A4 2n 0 0 
—2n O-A 0 0 
Mice i =a spen 
0 1 0 0-A 


The roots are A1,2= 0/434 = +2/n;the homogenous solution of system (9) is, 
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1 0 2n 2n 
0 0 —2in| -zint 2in | zint 
NH = Cy ol teh, tc i e + C3 e 
1 0 1 
1 0 2n 
E ; +t E E (cos 2nt — i sin 2nt) + c4 em (cos 2nt + i sin 2nt) 
1 0 1 
1 0 2n 0 
0 0 0 —2n]| . 
=c] jo +t 1 + (Cz + c3) o cos 2nt + | \sin 2nt 
1 0 1 0 
0 —2n 
+i (c2 — c3) E cos 2nt + : sin 2nt 
0 —1 
put c= (e+ c); B= I(&- c);then 
1 0 2n 0 0 
S 0 0 0 —2n| . —2n 
Xu = c 0 +t 1 + C2 0 cos 2nt + 1 sin 2nt ? + C3 1 cos 2nt + 
1 0 1 0 0 
—2n 
D 
o {sin 2nt (10) 
—1 


Where ci, c2, and cs are arbitrary linear independent constants. Then substitute into the 


system of Equations (2) we get, 


Xu = c iH +t lal} + Cz H cos 2nt + A sin 2nt} + C3 (t cos 2nt + El sin 2nt} 
(11) 


Now to obtain the particular- solutiony p, for the nonhomogeneous system (2) which is 


in the form, 


Xp = Ar(t) [$] +t [1 + A,(t) [$] cos 2nt + A sin ant} + A(t) (qu cos 2nt + [A] sin 2nt] (12) 


Since, the particular- solutionyp satisfied the system of Equations (2) then, 


asp da a 


where 2n i-a; 
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Xp = A(t) [i] +t D + AQ) (D + A,(t) [$] cos at + A sin at} + A(t) [$] (—a sinat) + A æ cos at} 


+ A3(t) H cos at + [A] sin at} + A;(t) (fl (—a sinat) + RA a cos at} 


= E o . fao [i +t pi +4A-(t) [$] cos at + M sin at) + A;(t) (fl cos at + LA] sin «ji 


= cfr] e eL] 4,0 {foe ae + [° 


E sin at} + A3(t) i] cos at + [5] sinc t} 


0 


Then 
A, (t) [$] +t [1 + Å, (t) [$] cos at + A sina t} + Å; (t) (qu cos at + ba sin at) +A,(t) (D 
t A2(t) lo] cos at + R 


-noge Lolla | 


| sin at) + A;(t) (L5.] cos at + [5 sina t} 
ee sin at} + A;(t) (L5.] cos at + L5] sina t} 


A,(t) [$] +t [1 + À;(t) [$] cos at + M sina t} + A;(t) (fl cos at + I] sin at} +A, = A,(t) lo] + 


0 1 
«51-9 
Then 
A, 
A, 
A; 
Equation (14) can be written as 

t A,(t) + sinat A,(t) + cosat A3(t) = (a — 1)A,(t) + Uy (15) 
Å (t) + cos at Å (t) — sinat A3(t) = —at A,(t) + Uy (16) 
where A, (t) is arbitrary function then, put 4, (t) = 0 and A, (t) = 0, this yields 
sin at Å (t) + cosat A3(t) = Uy (17) 
cos at A,(t) — sinat A3(t) = Uy (18) 
Multiply Equation (17) by sin at and Equation (18) by cos at and add then 


= A,() [77] + 


(14) 


1 cosat -—sinat 


f sinat  cosat 
0 0 0 














Ux 
Uy 
0 





(sin? at + cos? at) A,(t) = sinat Uy + cosat Uy 

A,(t) = sinat Uy + cosat Uy 

By integrating we get, 

A(t) = -Z cos at + "Y sin at (19) 
Multiply Equation (17) by cos at and Equation (18) by (— sin at )and add then 

(sin? at + cos? at) A3(t) = cosat Uy — sin at Uy 

A3(t) = cosat Uy — sinat Uy 

By integrating we get, 
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A3(t) = 7* sin at U, + Z cos at Uy (20) 
Substitute from Equations (19) and (20) into Equation (15) this yields 

tå (t) + sin at (Ecos at + Z sin at ) + cos at (sin at + cos at) = (a — 
1)A,(t) + Uy 

tÅ (t) — (a — 1)A,(t) = Uy —sinat (— cos at + 7 sin at) — cosat Cx sinat + 
Z cos at ) 

tÅ (t) — (a — 1)A,(t) = Uy — "x (— sin at cos at + sin at cos at ) — Z (sin? at + 
cos? at ) 


tå, (t) — 


(a-1) 
t 





A,(t) = (Ux + (21) 
Equation (21) represents linear ODE, which will be solved as follows 
d 1 1 Uy 
4; AiO i TE. = =z (Ux T P? ) 
By integrating to obtain the value of A1, 
1 U 1 
His cU Eee aep) (22) 


Now, substitute from Equations (19), (20) and (22) into Equation (12) the particular- 


solution is obtained 


Xp = (Ux — 5) (- —) iH +t mi T (— "X cos at + “sin at ) IH cos at + 





A sin at} + c sinat + x cos at ) (d cos at + E sin at} (23) 

Then the general analytical solutions for the nonhomogeneous system of ODEs (2) is given 
by 

X(t) = Xu + Xp (24) 


x(t) = E = c R +t [al} + C2 [i] cosat + A sin at} + C3 iH cos at + [A] sin at +} + (Ux — 


22 (- =) ([2] +t hi + (— "Ecos at c Z sin at ) (I2] cos at + M sin at) + Cx sinat + 


2 cos at ) (4 cos at + LA] sin at} (25) 





To apply the Variation of Parameters Method on a dynamical system, the locations of the 


libration points, Jacobi constant, and the stability of motion will be illustrated as follows. 


Locations of the Libration Points 

Let X = Ý =Z=0, and X = Y = Z = 0, and using Equation (25) at t = 0, to evaluate 
the initial condition. We put X(a) = Y (a) = b,, and X(a) = Ý (a) = Dp. 
When a 0, then, b, > 107^, b, > 0. We obtain the values of ci, cand c 
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a?10744-Uy _ @107*+Uy 


Which are given by c, = a(a-1) '"2 ^ a(a-1) 


and c, — 


These values are used through the system of Equations (2), this yield 


sx Cun a-8) _ EIS u AL zi (26) 
1 2 2 
2y  ü-)Xa-B)Y uY 3AnY — 
n°Y Sagat gas gee 0 (27) 


by solving Equation (26) numerically the locations of the collinear libration points are 
obtained, while the locations of the triangular libration points are calculated by solving 
Equations (26) and (27) together to obtain (X, Y) for each triangular point. 
The Jacobi constant 

It is well known that the Jacobi constant is given by Szebehely[19] 
C = (X? + Y?) - 2U (28) 
Since the Jacobi constant can be obtained at each libration point, then this enables to obtain 
the zero-velocity curves about each point. Since the variation of parameter solution of 
Equation (25) gives expressions for the (X, and Y) depending on time, which contains 
trigonometric functions depend on the angular velocity, these terms represent the short 
periodic orbits around the libration point understudy, the eccentricity and period of the 


orbit can be obtained by 





2-1 
e?-77 (29) 
21 
[E (30) 
2. 
where, c — IL S= Coefficient ofthe imaginary part of A Ibrahim[20]. 
i "XY 


Results and Discussion 

To apply the variation of parameter solution on the dynamical systems, the Sun-Earth- 
spacecraft system is considered. A Mathematica Cod is constructed to solve Equations (26) 
and (27) using Equation (25) to obtain the libration points for this model and to study the 
motion about each point. Then Table 1 illustrates the locations of libration points, Jacobi 
constant, eccentricity, and period of the orbit which are obtained also. 

The motion about each of the libration points is represented by showing the behavior of 
the phase space about the libration point understudy, the orbit obtained represents the 
short periodic orbits and it is obtained from Equation (25) with c1 = 0. Zero velocity curves 
show the regions at which the infinitesimal body moves with different values of energy 
levels. The first level is at the initial values of the Jacobi constant. Finally, the Poincare 


surfaces of section is a technique which represented the stability of the motion about the 
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libration points, it is a projection of the orbits from (X,Y,X) plane to (X,X) plane, each 
point represents an orbit about the libration point. 

Now, L1and L2 are chosen as an example of the collinear libration points, while L4 and L6 
are chosen as an example of the triangular libration points. 

From the results obtained Fig 2 shows periodic orbit about L1, Fig. 3 illustrates the zero- 
velocity curve and Fig.4 shows a Poincare surface about L1 and Fig 5 shows periodic orbit 
about L2, Fig. 6 illustrates the zero-velocity curve and Fig.7 shows a Poincare surface about 
L2. The same for L4, Fig 8 shows the periodic orbit about L4, Fig. 9 shows the zero-velocity 
curves and Fig.10 Poincare surface about L4. The same for L6 is illustrated in figures 11, 12, 
and 13. 


Table 1: The collinear and non-collinear libration points for the Sun-Earth system and their 
Jacobi constant C, eccentricity and period of the orbit about each point. 









































libration points Position Jacobiconstant  ccentricity | period ofthe orbit 
X Y 
L1 0.986675 0 2.98435 0.978755 3.13455 
L2 1.01201 0 2.98488 0.991976 3.13455 
L3 -1.00332 0 2.98463 1.6196 3.13455 
L4 0.986675 0.997497 3.3886 0.856303 3.13455 
L5 0.986675 -1.0025 3.39508 0.451298 3.13455 
L6 0.986675 9.00001*10^-10 2.98435 0.978755 3.13455 
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Fig.2: the phase space about L1 Fig. 3. Zero Velocity Curve at L1 , C=2.98435 


10 M.N. ISMAIL, A.H. IBRAHIM, SAHAR H. YOUNIS, GHADA F. MOHAMDIEN 














06 08 1.0 1.2 1.4 


Fig. 4. Poincare Surface at L1 
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Fig.5: the phase space about L2 Fig. 6. Zero Velocity Curve at L2 , C=2.98488 
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Fig. 7. - Poincare Surface at L2 
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Fig. 10: Poincare Surface at L4 
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Fig. 11: Phase Space at L6 Fig.12:Zero Velocity Curve at L6 ,C = 2.98435 
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Fig. 13: Poincare Surface at L6 
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Conclusion 

In this study, the Variation of Parameters Method is used to obtain the analytical 
solution of the restricted three-body problem. The obtained solution gives explicit 
expressions in X, Y depends on the time. The application of these solutions enables us to 
obtain the libration points and to study the stability of motion about these libration points, 
the results obtained are in a good agreement with the results obtained by Simmons [13], 
Kunitsyn[14], Schuerman[15] and Ibrahim!19 Therefore the variation of the parameter 


method is a good technique to be applied on the astro-dynamical systems. 
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